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(57) ABSTRACT 

Disclosed herein are systems, methods, and non-transitory 
computer-readable storage media for simulating propagation 
of an electromagnetic field, performing phase retrieval, or 
sampling a band-limited function. A system practicing the 
method generates transformed data using a discrete Fourier 
transform which samples a band-limited function f(x) with- 
out interpolating or modifying received data associated with 
the function f(x), wherein an interval between repeated copies 
in a periodic extension of the function f(x) obtained from the 
discrete Fourier transform is associated with a sampling ratio 
Q, defined as a ratio of a sampling frequency to a band-limited 
frequency, and wherein Q is assigned a value between 1 and 
2 such that substantially no aliasing occurs in the transformed 
data, and retrieves a phase in the received data based on the 
transformed data, wherein the phase is used as feedback to an 
optical system. 
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METHOD FOR UTILIZING PROPERTIES OF 
THE SINC(X) FUNCTION FOR PHASE 
RETRIEVAL ON 

NYQUIST-UNDER-SAMPLED DATA 

5 

BACKGROUND 

1 . Technical Field 

The present disclosure relates to phase retrieval and more 
specifically to sampling and reconstruction of the sinc(x) to 
function using phase retrieval. 

2. Introduction 

Broadly speaking, phase retrieval is a process used to 
retrieve an optical pupil phase and pupil amplitude based on 
images of a known obj ect received via an optical system. The 15 
phase and amplitude of an optical system are synonymous 
with aberrations, misalignments, or imperfections in the opti- 
cal system. Two general algorithm approaches are commonly 
utilized for phase retrieval. One approach is parametric based 
and the other is iterative-transform based. A number of varia- 20 
tions of both approaches have been developed, such as incor- 
porating diversity functions and one or more diversity 
images. However, these approaches introduce an aliasing 
effect in the resulting image data, and otherwise suffer from 
phase-wrapping discontinuities, ambiguous convergence to 25 
solutions, and estimation bias due to imperfect knowledge of 
the diversity function. What is needed in the art is an 
improved way to perform phase retrieval that avoids aliasing. 

SUMMARY 30 

Additional features and advantages of the disclosure will 
be set forth in the description which follows, and in part will 
be obvious from the description, or can be learned by practice 
of the herein disclosed principles. The features and advan- 35 
tages of the disclosure can be realized and obtained by means 
of the instruments and combinations particularly pointed out 
in the appended claims. These and other features of the dis- 
closure will become more fully apparent from the following 
description and appended claims, or can be learned by the 40 
practice of the principles set forth herein. 

Phase retrieval is a process used to recover an optical phase 
and amplitude, synonymous with aberrations, misalign- 
ments, or imperfections in an optical system, using images of 
a known object. The application of this technology can made 45 
to optical systems metrology, telescope control and align- 
ment, optical wavefront sensing and control, and can be per- 
tinent to adaptive optical systems for astronomical observing, 
security, and surveillance imaging. 

The approach disclosed herein takes advantage of sam- 50 
pling considerations for a band- limited function. The Fourier 
transform of this sampled, band-limited function is con- 
structed by periodic extension, i.e. spacing the copies in a 
definite way, such that minimal aliasing occurs for 1<Q<2. Q 
is the sampling ratio, which in turn is the ratio of the sampling 55 
frequency to the band-limited frequency. In optical systems, 
the sampling ratio is referred to as the image sampling param- 
eter and can be specified by Q=(lambda*f/#)/dx. Lambda is 
the assumed monochromatic wavelength of the scalar elec- 
tromagnetic field used to form the image being sampled, f/# is 60 
the f-number of the optical system, and dx is the image-plane 
sampling interval, or pixel size. 

There is a distinction in optics that electric fields should be 
sampled with Q>=1 to avoid aliasing, but that the irradiance 
measured by light detectors should be sampled with Q>=2 to 65 
avoid aliasing. This disclosure shows how phase retrieval can 
be performed with minimal aliasing for Q>=1 , despite the fact 


2 

that it makes use of irradiance data that would need Q>=2 to 
be interpolated without aliasing. 

The analysis considered here begins with the band-limited 
sinc(x) function and demonstrates that the sinc(x) function 
can be interpolated exactly using the Whittaker-Shannon 
sampling theorem for sampling ratios Q>=1. This interpola- 
tion is possible because of the extra space that exists between 
repeated copies of the sinc(x) function’s Fourier transform 
created by periodic extension. Using the band-limited prop- 
erty of the sinc(x) function and of the complex amplitude of 
the optical point-spread function, propagation of the electro- 
magnetic field can be simulated with no aliasing using the 
discrete Fourier transform (DFT) with 1<Q<2. Thus, phase 
retrieval can be performed with minimal aliasing on under- 
sampled point spread function (PSF) data, for sampling ratios 
2>Q>=1 . 

Disclosed are systems, methods, and non-transitory com- 
puter-readable storage media for simulating propagation of a 
monochromatic, scalar electromagnetic field, performing 
phase retrieval, and sampling a band-limited function. A sys- 
tem practicing the method generates transformed data using a 
discrete Fourier transform which samples a band-limited 
function f(x) without interpolating or modifying received 
data associated with the function f(x), wherein an interval 
between repeated copies in a periodic extension of the func- 
tion f(x) obtained from the discrete Fourier transform is asso- 
ciated with a sampling ratio Q, defined as a ratio of a sampling 
frequency to a band -limited frequency, and wherein Q has a 
value between 1 and 2 such that substantially no aliasing 
occurs in the transformed data, and retrieves a phase from the 
received data based on the received data and the transformed 
data. The phase can optionally be used as feedback to an 
optical system. The interval between repeated copies in a 
periodic extension of f(x) is associated with a sampling ratio 
Q, a ratio of sampling frequency to band-limited frequency. 
The system assigns Q a value between 1 and 2 so the discrete 
Fourier transform has substantially no aliasing, and retrieves 
a phase based on the received data and the transformed data. 
The phase can serve as feedback to tune, align, or otherwise 
correct an optical system. 

BRIEF DESCRIPTION OF THE DRAWINGS 

In order to describe the manner in which the above-recited 
and other advantages and features of the disclosure can be 
obtained, a more particular description of the principles 
briefly described above will be rendered by reference to spe- 
cific embodiments thereof which are illustrated in the 
appended drawings. Understanding that these drawings 
depict only exemplary embodiments of the disclosure and are 
not therefore to be considered to be limiting of its scope, the 
principles herein are described and explained with additional 
specificity and detail through the use of the accompanying 
drawings in which: 

FIG. 1 illustrates an example system embodiment; 

FIGS. 2A and 2B illustrate exemplary functions for con- 
tinuous sinc(v^x) and rect(v,v 6 )/v 6 for v b =2; 

FIG. 3 illustrates an exemplary continuous Fourier trans- 
form for sinc(v^x) for v b =2; 

FIGS. 4A and 4B illustrate an exemplary periodic exten- 
sion for Nyquist sampling for N=9 terms and N=5 1 terms; 

FIGS. 5A and 5B illustrate an exemplary periodic exten- 
sion for Nyquist sampling FFT of sampled sincfVpc); 

FIGS. 6A and 6B illustrate an example idealized periodic 
extension F 6p (v); 

FIG. 7 illustrates an example periodic extension f ipW; 

FIG. 8 illustrates an exemplary Whittaker- Shannon- 
Kotelnikov (WSK) interpolation for Q=l.l and N=23; 

FIG. 9 illustrates an exemplary weighted basis function for 
Q=l.l and N=23; 
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FIG. 10 illustrates a number of example aberration basis 
functions; 

FIGS. 11 A and 11B illustrate exemplary pupil, phase, and 
irradiance data for a 1-d simulation; 

FIG. 12 illustrates a first example method embodiment; 
and 

FIG. 13 illustrates a second example method embodiment. 

DETAILED DESCRIPTION 

Various embodiments of the disclosure are discussed in 
detail below. While specific implementations are discussed, it 
should be understood that this is done for illustration pur- 
poses only. A person skilled in the relevant art will recognize 
that other components and configurations may be used with- 
out parting from the spirit and scope of the disclosure. 

The present disclosure addresses the need in the art for 
performing phase retrieval without introducing aliasing 
effects into the results. A system, method and non-transitory 
computer-readable media are disclosed which simulate 
propagation of an electromagnetic field, perform phase 
retrieval, or sample a band-limited function. A system prac- 
ticing the method generates transformed data using a discrete 
Fourier transform which samples f(x) without interpolating 
or modifying received data associated with a band-limited 
function f(x). The interval between repeated copies in a peri- 
odic extension of f(x) is associated with a sampling ratio Q, a 
ratio of sampling frequency to band-limited frequency. A 
discussion of a basic general purpose system or computing 
device in FIG. 1 which can be employed to practice the 
concepts is disclosed herein. A more detailed description of 
the methods and algorithms, and other variations will then 
follow. The disclosure now turns to FIG. 1. 

With reference to FIG. 1, an exemplary system 100 
includes a general-purpose computing device 100, including 
a processing unit (CPU or processor) 120 and a system bus 
110 that couples various system components including the 
system memory 130 such as read only memory (ROM) 140 
and random access memory (RAM) 150 to the processor 120. 
The system 100 can include a cache 122 of high speed 
memory connected directly with, in close proximity to, or 
integrated as part of the processor 120. The system 100 copies 
data from the memory 130 and/or the storage device 160 to 
the cache 122 for quick access by the processor 120. In this 
way, the cache 122 provides a performance boost that avoids 
processor 120 delays while waiting for data. These and other 
modules can be configured to control the processor 120 to 
perform various actions. Other system memory 130 may be 
available for use as well. The memory 130 can include mul- 
tiple different types of memory with different performance 
characteristics. It can be appreciated that the disclosure may 
operate on a computing device 100 with more than one pro- 
cessor 120 or on a group or cluster of computing devices 
networked together to provide greater processing capability. 
The processor 120 can include any general purpose processor 
and a hardware module or software module, such as module 
1 162, module 2 164, and module 3 166 stored in storage 
device 160, configured to control the processor 120 as well as 
a special-purpose processor where software instructions are 
incorporated into the actual processor design. The processor 
120 may essentially be a completely self-contained comput- 
ing system, containing multiple cores or processors, a bus, 
memory controller, cache, etc. 

The system bus 110 may be any of several types of bus 
structures including a memory bus or memory controller, a 
peripheral bus, and a local bus using any of a variety of bus 
architectures. A basic input/output (BIOS) stored in ROM 
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140 or the like, may provide the basic routine that helps to 
transfer information between elements within the computing 
device 100, such as during start-up. The computing device 
100 further includes storage devices 160 such as a hard disk 
5 drive, a magnetic disk drive, an optical disk drive, tape drive 
or the like. The storage device 160 can include software 
modules 162, 164, 166 for controlling the processor 120. 
Other hardware or software modules are contemplated. The 
storage device 160 is connected to the system bus 110 by a 
l drive interface. The drives and the associated computer read- 
able storage media provide nonvolatile storage of computer 
readable instructions, data structures, program modules and 
other data for the computing device 100. In one aspect, a 
hardware module that performs a particular function includes 
the software component stored in a non-transitory computer- 
1 5 readable medium in connection with the necessary hardware 
components, such as the processor 120, bus 110, display 170, 
and so forth, to carry out the function. The basic components 
are known to those of skill in the art and appropriate variations 
are contemplated depending on the type of device, such as 
20 whether the device 100 is a small, handheld computing 
device, a desktop computer, or a computer server. 

Although the exemplary embodiment described herein 
employs the hard disk 160, it should be appreciated by those 
skilled in the art that other types of computer readable media 
25 which can store data that are accessible by a computer, such as 
magnetic cassettes, flash memory cards, digital versatile 
disks, cartridges, random access memories (RAMs) 150, read 
only memory (ROM) 140, a cable or wireless signal contain- 
ing a bit stream and the like, may also be used in the exem- 
30 plary operating environment. Non-transitory computer-read- 
able storage media expressly exclude media such as energy, 
carrier signals, electromagnetic waves, and signals per se. 

To enable user interaction with the computing device 100, 
an input device 190 represents any number of input mecha- 
35 nisms, such as a microphone for speech, a touch-sensitive 
screen for gesture or graphical input, keyboard, mouse, 
motion input, speech and so forth. An output device 170 can 
also be one or more of a number of output mechanisms known 
to those of skill in the art. In some instances, multimodal 
40 systems enable a user to provide multiple types of input to 
communicate with the computing device 100. The commu- 
nications interface 180 generally governs and manages the 
user input and system output. There is no restriction on oper- 
ating on any particular hardware arrangement and therefore 
45 the basic features here may easily be substituted for improved 
hardware or firmware arrangements as they are developed. 

For clarity of explanation, the illustrative system embodi- 
ment is presented as including individual functional blocks 
including functional blocks labeled as a “processor” or pro- 
50 cessor 120. The functions these blocks represent may be 
provided through the use of either shared or dedicated hard- 
ware, including, but not limited to, hardware capable of 
executing software and hardware, such as a processor 120, 
that is purpose-built to operate as an equivalent to software 
55 executing on a general purpose processor. For example the 
functions of one or more processors presented in FIG. 1 may 
be provided by a single shared processor or multiple proces- 
sors. (Use of the term “processor” should not be construed to 
refer exclusively to hardware capable of executing software.) 
60 Illustrative embodiments may include microproces sor and/ or 
digital signal processor (DSP) hardware, read-only memory 
(ROM) 140 for storing software performing the operations 
discussed below, and random access memory (RAM) 150 for 
storing results. Very large scale integration (VLSI) hardware 
65 embodiments, as well as custom VLSI circuitry in combina- 
tion with a general purpose DSP circuit, may also be pro- 
vided. 
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The logical operations of the various embodiments are 
implemented as: (1) a sequence of computer implemented 
steps, operations, or procedures running on a programmable 
circuit within a general use computer, (2) a sequence of com- 
puter implemented steps, operations, or procedures running 5 
on a specific-use programmable circuit; and/or (3) intercon- 
nected machine modules or program engines within the pro- 
grammable circuits. The system 100 shown in FIG. 1 can 
practice all or part of the recited methods, can be a part of the 
recited systems, and/or can operate according to instructions 10 
in the recited non-transitory computer-readable storage 
media. Such logical operations can be implemented as mod- 
ules configured to control the processor 120 to perform par- 
ticular functions according to the programming of the mod- 
ule. For example, FIG. 1 illustrates three modules Modi 162, 15 
Mod2 164 and Mod3 166 which are modules configured to 
control the processor 120. These modules may be stored on 
the storage device 160 and loaded into RAM 150 or memory 
130 at runtime or may be stored as would be known in the art 
in other computer-readable memory locations. 20 

Having disclosed a basic computing device which can 
practice the method, the disclosure now turns to a discussion 
of sampling and reconstruction of the sinc(x) function. 

For the propagation of electromagnetic fields calculated as 
the Fourier transform of a general pupil geometry, this 25 
approach can exploit some special qualities of the electric 
field Fourier transform, F, to minimize aliasing, or non- 
uniqueness of the Fourier transform. In particular, the Fourier 
transform of a sampled, band-limited function can be made 
periodic by applying a “periodic extension” in Fourier analy- 30 
sis. Thus the sampled, band-limited function can be expressed 
as a Fourier series. The periodic replicas of F have minimal 
overlap, and therefore exhibit no aliasing for values of Q>=1 . 
This result differs from the conventional viewpoint in phase 
retrieval that Q must be greater than or equal to 2 when the 35 
data (images, which are measurements of the irradiance at the 
detector) are not modified or interpolated. This disclosure 
demonstrates how the Fourier transform of a sampled, band- 
limited optical function can be periodically extended in a way 
that minimizes aliasing. Thus, phase retrieval can be imple- 40 
mented in a way that minimizes aliasing by an appropriate 
choice of sampling variables for constructing the Fourier 
transform by periodic extension. 

The sampling theorem can be derived in a way that empha- 
sizes two assumptions of the theorem explicitly. For example, 45 
the theorem can be expressed in terms of two length scales 
that are derived from the data sampling frequency, v A , and the 
data band-limited frequency v b . With the substitution Q=v^/ 
v A , the results can be expressed by the following function: 


= £ ^ f( X n)^c[2 Vb (x n ~X)] 
n=-oo 

or 

+oo 

f(x)\ Q ^ 2 = 2 f(x n )sinc(n-xl&x) 

n=-oo 

Too 

= 2 fMsmc[v A (x n -x)] 


Equation (2), with Q=2, is the standard form of Equation 
(1) (Whittaker-Shannon sampling theorem). The derivation 


50 

(1) 

55 

(2) 

60 

65 
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leading to Equation (1) can lead to an alternative way of 
looking at the problem. The preferred basis set for interpola- 
tion can be found by varying the frequency component of the 
basis functions to minimize the side lobe artifacts in the 
superposition of the weighted basis functions. Some 
examples comparing Equations (1) and (2) for the “canoni- 
cal” sinc2 (x) function can illustrate the results. However, 
special functions exist where QE(1,2] is perfectly valid and 
no aliasing occurs. Indeed, one such example is the sinc(x) 
function. Examining the sinc(x) function is instructive 
because its coverage in frequency space is half that of the 
sinc 2 (x) function, and therefore, some subtleties exist in the 
relationship between its sampled bandpass and its continuous 
bandpass. By making this distinction, the generality of the 
analysis is made explicit. Some numerical examples are also 
discussed to give some basic insight into the sampled sinc(x) 
function in the context of the Whittaker- Shannon-Kotelnikov 
(WSK) assumptions. 

Some further examples of “sampling” in this context are of 
interest from the perspective of scalar diffraction theory since 
the focal plane electric field can be modeled in each of the two 
transverse dimensions using the sinc(x) function. 

The “continuous” results emphasize the role of sampling. 
The convention adopted here is that x is the spatial variable 
and the sinc(v^,x) function is defined by 


fix) = sinc(vfc*) = 


sin(7TV£*) 

7lV b X 


(3) 


where \ b is the band-limited frequency specifying the non- 
zero extent of the Fourier transform. The Fourier transform 3 
of Equation (3) is given by the scaled rect(x) function below: 


Ftiv) = Jismcivtx)} (4) 

1 

= — rect(v , v b ) 

Vb 

_ ( l/v t , for | v | < v b j2 
( 0, otherwise 


with plots of Equations (3) and (4) shown in FIGS. 2 A and 2B 
for v b =2, noting that the spatial period associated with \ b in 
this case is given by 2/v 6 =l. FIGS. 2A and 2B illustrate 
exemplary functions for continuous sinc(v^x) and rect(v,v 6 )/ 
\ b for v b =2. A slightly different notation is used here for the 
rect function than the standard form. Specifically, a two- 
parameter form, rect(x, y), is used over the usual notation of 
rect(x/y) to avoid confusion when scaling and using dimen- 
sionless variables. 

For reference, the Fourier transform convention adopted 
here is shown below: 


F(v) = . 3 {/(*)} 


■£ 


dxf(x)e 1 


(5) 


and its inverse: 


/W = 3- 1 {F(v)} = 


£ 


dvF(vy 2 ™ 


( 6 ) 
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The main assumptions from classical sampling theory are 
summarized succinctly in terms of two main conditions. First, 
f(x) is band-limited: 

F 6 (v)=3(/(x)}= 0 for v<£[-v b ,v b ] and lv 6 l<co (7) 5 

Second, F 6 (v) is constructed by periodic extension over 
the sampling interval v A : 

F bp (v)=F b {v+nv A ) , with v a =1/Ax, for VvGtK. and n=0, 

1 , 2 ,... ( 8 ) 

1C 

By application of Equation ( 8 ), F^(v) can be calculated 
using a Fourier series expansion in v A as shown below: 


+°° +°° (9) i « 

Fbp(v) = 2 e„e-' 2 ™ (1/v a ) = ^ 

n=—oo n=—oo 


The construction in Equation (9) is not computable in 
practice and so when using a finite number of terms the 20 
discrete Fourier series (DFS) below 


F<wv(v) = 2 c„c- iW( ” )/v A s Yj c » e ' 


ilnvn' (n)Ax 


( 10 ) 


where 


n'{n)=\--{N -l). 


’ + 2 (N " 1} 


(ID 


30 


maps the index range over N total points including values of 
xOforxi^O, n'(n) can be defined as n'(n)=[0, 1,2,... ,N-1]. 
Note that Equation (1 1 ) is not the only choice. Another option 
is given below: 35 


n'(n)= -^N, 


■ + 2 (A, - 1) 


( 12 ) 


The disclosure now turns to a periodic extension of the 
sampled sine function. By inspection of the exemplary con- 40 
tinuous Fourier transform for f(x)=sinc(v^x) for v^=2 as 
shown in FIG. 3, F^(v)=rect(v, v b )/v b is only non-zero over 
W b \^\ b l2 and is the continuous result. In sampling the f(x) 
function, v A is introduced which determines how the repeated 
copies of the F^(v) are placed with respect to v b . So when 45 
constructing the F 6 (v) by periodic extension, or when f(x) is 
physically sampled (such as when f(x) is a measured data 
signal and is sampled by a physical detection process) and the 
Fourier transform is then applied to this sampled data, the 
interval for periodic extension is defined by v A =Q v b . This 
interval defines the spacing between the repeated copies of 50 
the F 6 (v), and is helpful to understanding the true bandpass, 
v b , of the sinc(v 6 x) function, and its apparent bandpass, v' b . 
Equation (10) can be alternatively expressed in terms of Q: 


55 

N - 1 N-l (13) 

Fipjv(v) = Yj ^ c n e-^ W(n)x b/Q 

n = 0 n = 0 


Continuing the calculation produces the Fourier coeffi- 60 
cients c„ in Equation (10) that are given by (see Appendix A): 


sin c[njQ) 
va 


(14) 

65 
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and therefore the “periodically extended” rect(v, v b )/v b func- 
tion is approximated for finite N as 


1 


F bp( v ) = — + 
V A 



sinefn / Q)cos(2jmv j v a ) 


(15) 


To illustrate, consider the situation for Nyquist sampling 
and arbitrarily choose v b =2: 

Q=v A /v h =v A /2=2 => v A =4 Ax = Va (16) 

FIGS. 4A and 4B illustrate an exemplary periodic exten- 
sion for Nyquist sampling for N=9 terms and N=51 terms. 
The DFS expansion shown in Equation (1 5) is shown in FIG. 
4 A for N=9 terms, and for comparison, the approximation 
shown in FIG. 4B is for N=51 terms. A prominent feature is 
the oscillating nature of the Fourier series due to the trunca- 
tion of terms in going from Equation (9) to Equation (10). 
This ringing effect is the Gibbs phenomena. In summary of 
FIGS. 4Aand 4B, the periodic extension of F 6 (v) is calculated 
using the DFS via Equation (10). 

Equivalently, the transform of the data is calculated 
directly by Fourier transforming the sampled sinc(v 6 , x) func- 
tion using the Fast Fourier Transform (FFT), rather than by 
the DFS expansion. The DFT can also be used to derive from 
the fundamental period of the DFS. For comparison, this 
calculation is shown in FIGS. 5A and 5B. FIGS. 5 A and 5B 
illustrate an exemplary periodic extension for Nyquist sam- 
pling FFT of sampled sinc(v^x), noting that the Gibbs ringing 
can be made arbitrarily small (but never zero) as discussed 
below with respect to minimal aliasing for Q« 1 . 

Next the disclosure turns to a discussion of special func- 
tions that exist where QE(1,2] which are perfectly valid and 
in which no aliasing occurs. FIGS. 6 A and 6 B illustrate an 
example idealized periodic extension F^, (v). FIG. 6 A shows 
a Nyquist sampled version of the idealized periodic extension 
F^(v). This idealized construction shows that no overlap 
occurs between neighboring copies of the repeated F^(v), for 
all Q>=1 , therefore no aliasing is present. Note, however, that 
this statement is only valid for the idealized case with no 
truncation error. In practice, the data is “windowed” or trun- 
cated leading to Gibbs “ringing” and therefore no matter how 
small, these errors always “interfere” at some level in the 
space between the gaps of the F^(v). To further illustrate 
Equation (10) (or its specific result, Equation (15)), the cal- 
culation is shown in FIG. 7 for Q=1 .1 and N=1 1 . 

Q can be made as close to 1 as desired, limited only by 
available computer memory or other computational factors. 
Q in turn dictates the window size of the data. Truncation can 
also occur due to the finite extent of the detector. Some 
additional implications are discussed below. 

The disclosure now turns to Whittaker-Shannon-Kotelni- 
kov (WSK) reconstruction. As disclosed above, the sampled 
sinc(v^x) function can be interpolated to arbitrary precision 
for Q>=1 . Nevertheless, it is instructive to demonstrate the 
reconstruction for a specific numerical example. Continuing 
with the example choosing Q=l.l and v 6 =2, the associated 
quantities for the data sampling frequency and the sample 
spacing are 


Q=v A /v b =v A /2=l. 1 =>v a =2.2=>A*~0.45 (17) 

The WSK interpolation formula shown in Equation (2) is 
then applied to reproduce f(x) from a set of N=23 uniformly 
sampled data values. The sample points of f(x) are x„=n Ax 
and marked by the vertical dashed lines in FIG. 8 . The WSK 
result produces the interpolated function values, j I (x , m ), 
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shown in the figure at points x , w =mAx' for 221 points for 
lx I <=2.5, marked by circles and separated along the x-axis. 

The weighted sinc(x) basis functions are distributed across 
a number of the data sample points and are shown in FIG. 9. 
The sum of the data sample points in FIG. 9 gives the result 
shown in FIG. 8. In summary, a Q>=1 sampled sinc(v 6 x) 
function can be reconstructed to arbitrary precision. FIG. 8 
shows that a good approximation can be achieved using just a 
modest number of 23 data sample points. The RMS error 
between f(x)=sinc(v 6 x) and its interpolated counterpart is 
calculated in units normalized to one, as shown below: 


D=1 . In order to derive this particular basis set, the system can 
start with the set of “seed” polynomials, which for simplicity 
in this example are non-negative integer powers of the spatial 
frequency variable, v". Then the system generates higher 
5 spatial frequency polynomials using the Gram-Schmidt 
orthogonalization procedure over the interval 


to 



li 

+ 2 . ' 


RMS[f(x)-f^)]=2.6xlO- 3 (18) 

The disclosure will now provide a phase-retrieval example. 
The sinc(v^x) function has be interpolated to arbitrary accu- 15 
racy using the WSK theorem for all Q>=1 . FIGS. 6 A and 6B 
shows that this ability is due to the extra “space” that exists 
between copies of the neighboring, F 6 (v)=rect(v,v 6 )/v 6 , cre- 
ated by periodic extension, which in turn is implied by sam- 
pling. Similarly, this observation suggests that phase retrieval 20 
can be performed on under-sampled data, with 1<=Q<2, by 
identifying the focal plane electric field in each transverse 
dimension with sinc(v 6 x) function. The disclosure will dem- 
onstrate some of these details and results for an iterative- 
transform class of phase-retrieval algorithms. 25 

The application to phase retrieval proceeds by identifica- 
tion of sinc(v^x) with the focal plane electric field in each 
transverse dimension. Thus, the Fourier conjugate of the exit 
pupil amplitude is F 6 (v)=rect(v,v 6 )/v 6 in each transverse 
dimension. In one implementation, the phase-retrieval algo- 30 
rithm uses a discrete Fourier transform (DFT) propagator in 
combination with an iterative-transform type phase-retrieval 
algorithm, specifically an implementation of the Misell-Ger- 
chberg- Saxton algorithm. A phase retrieval approach using 
values of Q>«1 can be used. 35 

The discrete or “aliased” version of the Fourier transform 
can be derived from the sampling theorem as the fundamental 
period of the DFS assuming that v A »l, giving F^(v) rather 
than the F^ (v). The transform can be defined at M interpo- 
lation points by a sum over N data samples: 40 


n- i (19) 

F b (v m ) * DFT(v m ) = /(*„)e-' 2 ’ r '" WW 

« =1 45 

where \ m denotes sample frequencies over the fundamental 
period of the F*p(v): 

50 


v m e 



Va 




( 20 ) 


The resulting basis set is listed in Appendix B, noting that the 
analog of these one dimensional polynomials to their two 
dimensional Zemike counterparts are listed in FIG. 10 and in 
Appendix B. 

These v w seed polynomials yield the Legendre polynomials 
when applying the Gram-Schmidt orthogonalization over the 
interval vE[-l, 1], rather than over 


v e 



n 

+ 2 . ’ 


The polynomials defined in Appendix B are not identical to 
the Legendre polynomials but are closely related since the 
change of variable, v-*v ,= 2v, leads to the familiar form over 
the interval vE[-l, 1]. Exemplary code for generating the 
polynomials using Wolfram Research’s Mathematica™ soft- 
ware is given in Appendix C. The one-dimensional aperture 
function can be modeled using the rect function: F 6 (v)=rect 
(v,v 6 ), with v b =1 /2 and 


is a dimensionless spatial frequency variable and D is the 
aperture diameter. 

The disclosure now turns to a discussion of the data for 
Q>«1. Given the basis functions defined in Appendix B and 
the observations above regarding minimal aliasing for the 
sinc(v^x) function, phase retrieval performance can be dem- 
onstrated with negligible aliasing for an arbitrary value of 
Q>«1. For example, let Q=1.06. The exact value, Q=l, is 
problematic because the replicated copies of the F 6 (v) are 
coincident at this “asymptotic” numerical value, as shown in 
FIG. 6A. Thus, an interpretation of aliasing is difficult at this 
value, so Q is confined to values greater than 1 for the present 
discussion for clarity. Some phase-retrieval approaches use 
values of Q<1 . The system can generate a phase with the 
following algorithm: 


For reference, the inverse transform of Equation (1 9) is 55 


n- 1 (21) 

f(x m ) * DFT\F ) = (MAX)- 1 ^] F b (v n )e +i2KXm(nAv) 

n= 0 

60 

where x a =1/Av is the function spatial period, not to be con- 
fused with the function sampling period Ax. 

A one-dimensional phase basis set and aperture can facili- 
tate the visualization. The basis polynomials, L k (v), are 65 
defined as the set of polynomials orthogonal on the unit 
one-dimensional rectangular aperture, i.e., with diameter, 


10 (22) 

<p(v) = ^ C nU(v) 

n= 1 

by choosing the basis coefficients c„ randomly from a 
Gaussian distribution and the L„(v) correspond to the basis set 
in Appendix B. The basis coefficients, c", can be displayed for 
a single realization and can be used as a consistency check 
with the phase fitting procedure by applying a least-squares 
decomposition on this starting phase. 

The system can generate irradiance data using this phase 
realization for the diversity defocus value, c 3 =LX (where 
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1 X=2jt, noting that the particular value of c 3 can be chosen 
with some latitude, the goal being to have the image data fill 
a reasonable portion of the detector. From Appendix C, the 
quadratic diversity-defocus term is given by 


G .diversity = V^V 2 - ij 


(23) 


In practice the value for c 3 can be chosen to enhance phase- 
retrieval estimation performance based on the spatial fre- 
quency content of (|). Example one-dimensional pupil, phase, 
and irradiance data for Q=1 .06 are calculated and shown in 
FIGS. 11 A and 11B. To specify Q, the system can consider a 
Nyquist sampled point spread function (PSF) with Airy diam- 1 5 
eter=4.88 pixels across. When the Nyquist sampling is arbi- 
trarily set to N a ^=512 pixels, for Q«1.06 the under-sampled 
data set corresponds to N data pixels: 


A I hlyn 

Ndata = (2yfj = 272 pixels 


20 

(24) 


To be consistent with the earlier sine function examples, 
the variable x in sinc(v^x) applies to the spatial domain (focal 25 
plane) while the Fourier domain is labeled using a dimension- 
less spatial frequency variable in the pupil (as noted in FIG. 
10): F 6 (v)=rect(v,v 6 )/v 6 . The DFT in Equation (19) can then 
be used to calculate the irradiance using 272 sample points. 
The implied number of pupil samples (wavefront) are thus by 30 
N ., 7=256, resulting in the following: 


Ndata (25) 


Using the under-sampled irradiance data of FIGS. 11 A and 
11B as input to the phase-retrieval algorithm, the system 
executes 200 iterations and compares the output of the phase- 
retrieval result to the input phase. The true phase coefficients 40 
(simulated by the computer algorithm) can be compared to 
the recovered coefficients using a least-squares decomposi- 
tion. A small piston term L x (v) is recovered, which is char- 
acteristic of the iterative-transform approach. The piston term 
causes a small offset in the true phase compares to the recov- 45 
ered phase, which accounts for the majority of the overall 
wavefront difference between the two. The root mean square 
(RMS) difference between the actual and recovered phase 
values is 

RMS error=rms (<b^ -<IW J~8x 1 0~ 4 \ (26) 50 

and can be driven to even smaller values with additional 
iterations, yielding an error term dominated by piston. 

The example results discussed herein include no other 
noise or detector effects and serve mainly to illustrate that no 
significant limitation prevents the phase-retrieval algorithm 55 
from successfully recovering phase values from under- 
sampled data for Q>=1 . 

One significant advantage of using discrete Fourier trans- 
forms (DFTs) rather than fast Fourier transforms (FFTs) in 
phase retrieval is that the phase retrieval is performed directly 6Q 
at Q«1.06. In other words, the DFT calculates the under- 
sampled model values directly, without the need for further 
interpolation or modification of the data. This approach can 
be contrasted to the usual approach to handling under- 
sampled data in phase retrieval: the data is first interpolated to 
Nyquist sampling and then the FFTs are implement using 65 
minimum pad-sizes that are equal to 2 times the number of 
pupil samples. In this regard, the DFT provides a more flex- 


ible approach, which has the advantage of requiring no modi- 
fication of the original data values through interpolation. 
These various data-modifying procedures can introduce spu- 
rious and possibly non-physical artifacts into the retrieved 
phase. 

Having disclosed some basic system components, the dis- 
closure now turns to the exemplary method embodiments 
shown in FIGS . 12 and 13 . For the sake of clarity, the methods 
are discussed in terms of an exemplary system such as is 
shown in FIG. 1 configured to practice the methods. 

FIG. 12 illustrates a first example method embodiment for 
sampling a band-limited function. The system 100 generates 
transformed data using a discrete Fourier transform which 
samples a band-limited function f(x) without interpolating or 
modifying received data associated with the function f(x), 
wherein an interval between repeated copies in a periodic 
extension of the function f(x) obtained from the discrete 
Fourier transform is associated with a sampling ratio Q, 
defined as a ratio of a sampling frequency to a band-limited 
frequency, and wherein Q is assigned a value between 1 and 
2 such that substantially no aliasing occurs in the transformed 
data (1202). The band-limited function can be a sinc(x) func- 
tion. In one variation, Q=(lambda * f/ #)/ dx , such that lambda 
is a monochromatic wavelength of a scalar electromagnetic 
field, f/# is an f-number of an optical system, and dx is an 
image plane sampling interval. In one aspect, the received 
data was collected via a confined circular aperture of diameter 
D. 

The system 100 retrieves a phase in the received data based 
on the transformed data, wherein the phase is used as feed- 
back to an optical system (1204). The system 100 can further 
adjust the optical system based on the feedback or the optical 
system can auto-adjust based on the feedback. 

FIG. 13 illustrates an alternative second example method 
embodiment for sampling a band-limited function. The sys- 
tem 100 receives data such as from an aperture or point source 
(1302), samples the received data, generating a periodic 
expansion representation of the received data (1304), sets the 
Q value between 1 and 2 to avoid aliasing (1306), and 
retrieves the phase without interpolation and without modi- 
fication of the original received data (1308). 

Embodiments within the scope of the present disclosure 
may also include tangible and/or non-transitory computer- 
readable storage media for carrying or having computer-ex- 
ecutable instructions or data structures stored thereon. Such 
non-transitory computer-readable storage media can be any 
available media that can be accessed by a general purpose or 
special purpose computer, including the functional design of 
any special purpose processor as discussed above. By way of 
example, and not limitation, such non-transitory computer- 
readable media can include RAM, ROM, EEPROM, CD- 
ROM or other optical disk storage, magnetic disk storage or 
other magnetic storage devices, or any other medium which 
can be used to carry or store desired program code means in 
the form of computer-executable instructions, data structures, 
or processor chip design. When information is transferred or 
provided over a network or another communications connec- 
tion (either hardwired, wireless, or combination thereof) to a 
computer, the computer properly views the connection as a 
computer-readable medium. Thus, any such connection is 
properly termed a computer-readable medium. Combinations 
of the above should also be included within the scope of the 
computer-readable media. 

Computer-executable instructions include, for example, 
instructions and data which cause a general purpose com- 
puter, special purpose computer, or special purpose process- 
ing device to perform a certain function or group of functions . 
Computer-executable instructions also include program 
modules that are executed by computers in stand-alone or 
network environments. Generally, program modules include 
routines, programs, components, data structures, objects, and 
the functions inherent in the design of special-purpose pro- 
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cessors, etc. that perform particular tasks or implement par- 
ticular abstract data types. Computer-executable instructions, 
associated data structures, and program modules represent 
examples of the program code means for executing steps of 
the methods disclosed herein. The particular sequence of such 5 
executable instructions or associated data structures repre- 
sents examples of corresponding acts for implementing the 
functions described in such steps. 

Those of skill in the art will appreciate that other embodi- 
ments of the disclosure may be practiced in network comput- 
ing environments with many types of computer system con- 10 
figurations, including personal computers, hand-held 
devices, multi-processor systems, microprocessor-based or 
programmable consumer electronics, network PCs, mini- 
computers, mainframe computers, and the like. Embodi- 
ments may also be practiced in distributed computing envi- 15 
ronments where tasks are performed by local and remote 
processing devices that are linked (either by hardwired links, 
wireless links, or by a combination thereof) through a com- 
munications network. In a distributed computing environ- 
ment, program modules may be located in both local and 
remote memory storage devices. 20 

The various embodiments described above are provided by 
way of illustration only and should not be construed to limit 
the scope of the disclosure. For example, the phase retrieval 
principles disclosed herein can be applied in optical systems 
metrology, telescope control and alignment, and optical 25 
wavefront sensing and control, and is pertinent to adaptive 
optical systems for astronomical observing, security, and sur- 
veillance imaging. Those skilled in the art will readily recog- 
nize various modifications and changes that may be made to 
the principles described herein without following the 
example embodiments and applications illustrated and 
described herein, and without departing from the spirit and 
scope of the disclosure. 


APPENDIX A 

Periodic Extension by Fourier Series 


35 


Periodic Extension of the sampled sinc(v 6 x) transform is 
accomplished using the Fourier series expansion in Equation 
(9). Therefore, it is instructive to derive this result as it also 40 
helps to illustrate the aliasing errors introduced by the trun- 
cation of terms. To begin, multiply Equation (9) by e z2jtWvA 
and integrate both sides over one full F^(v) function period, 
defined from -v A / 2 to +v A /2. Interchanging the order of inte- 
gration and summation gives 45 


dvF^v)^^ = J\„jT 


»A/2 


(Al) 


dve i27:v{m ~ n)Ax 


-va/2 


50 
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then evaluating the right hand side of Equation (Al): 


+= (A2) 

Va j c«sin c(m-n) = v A c m 

n=—oo 


Combining Equation (A2) with the left hand side of Equa- 
tion (Al): 


1 r +v A/2 1 r+v a/2 (A3) 

c n =— dvFUv)e i2mv/v * = — dvF b ( V y 2nnv/v a 

VaJ_ V a /2 V A J_ Va/2 

using the fact that over a single function period, F^(v)=F^(v). 
But V b (y)=V 2 =\N b (from Equation (4)) is only nonzero over 
-v b !2 to +v 6 /2, which leads to Equation (15): 


c _ _J_ ^v/v A _ si nprrc/6) _ sin cjnlQ) (A4) 

" ” V A J-v a /2 v V b ) 7mv b V A 


The general Fourier series expansion is expressed in the usual 
form: 


a o v— i 

F bp (v ) = — + 2_j [a n cos(27mv/ v A ) + b n (sin(2jmv / v A )\ 


(A5) 


but because F^(v)=rect(v,v^)/v 6 is an even function, it follows 
that b„=0, a w =c„+c_ w =2c w and thus 


2 . 2 
: — sin c(n/ Q), with a® = — 
va v a 


(A6) 


Substituting these values into Equation (A5) yields the 
following equation: 


1 2 ^ 

F bp (v ) = — + — ) sinc(re/Q)cos(27rov/v A ) 
p v A v A ^ 


(A7) 


APPENDIX B 
Polynomial Basis Functions 


Term # Polynomial : L^( v ) : Zernike Analog 

1 1 Piston 

2 V3(2 v) Tilt 

3 2 Defocus 

77(6.--) 

4 V7(20 v 3 -3 v) Coma 

V9"(70v 4 - 15v 2 + ? j 


5 


Spherical 
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-continued 


Term # 

Polynomial: L t (v ): 

Zernike Analog 

6 

VTT^252v 5 -70v 3 + ^vj 

Astigmatism 

7 

VTT(924v 6 - 315v 4 + ^ v 2 - ^ j 

2 W order spherical 

8 

VTT(3432v 7 - 1386v 5 + ^v 3 -^v] 

2” rf order astigmatism 

9 

, — ( -8 -6 3465-4 315-2 35 \ 

VTT 12780v -6006v + — — v - — v + 

v 48 128 ) 

3 rd order spherical 

10 

, — f „ 9 -7 9009-s 1155-3 315 A 

VTT 48620v -25740v + — v - — v + — vl 

3 rd order astigmatism 


APPENDIX C 


20 

Calculate Basis: 


Wolfram Research’s Mathematica™ Code 

One Dimensional Basis Functions 
Orthonormalized Polynomials 


n_, v_, vl_, v2_] = 


*p[n, v] 


tf/[n, v ] 


Integrate^ [«, v]2 ’ (v ’ vl ’ v21] 


f 


«[*> v]tf/[j, v ] 


Intermediate Coefficients 

a[ i_, j_, v_] = -Integrate 

[«[/, v]<p[j, v, vl, v2], {v, vl, v2}] - 

Inner Product 


ip[i_, j-1 = Integrate^/, v, vl, v2]<p[j, v, vl, v2], {v, vl, v2}] 

ip\i, v]| (/[j, v] 


25 basis = {}; terms = 10; 

Do[ 

W[m,v] = u[m,v] + Sum[a[m,k,v] <|)[k,v,vl,v2],{k,0,m - 1}]; 
AppendTo[basis,Collect[(|)[m,v,vl ,v2] ,v]],{m,0, terms- 1 }] ; 
30 TableForm [basis] 


1 

2V3v 

35 

+ 6V3"v 2 


40 -3V7v + 20V7v 3 


■dv 8 


-45v 2 + 210v 4 


45 


15V1U 


- 70 VTT" v 3 + 252VTT" v 5 


I 


J ^ lAE v] 2 Jv^ JjVty, v] 2 rfv 


: V 


50 


55 


5 VTT 105 VTT v 2 


16 


-315VTTv 4 +924VTTv 6 


35 VTT v 315VTTV 


- 1 386 VTT V 5 + 3432 VTT v 7 


35VTT 315VT7 V 2 3465 VTT V 4 , — , , — s 

— — 5 + 6006VTT v 6 + 12870VTT v 8 

128 8 4 


Starting Basis (Assume Weighting Function, W(v)=l) 60 

u[n_,v_ ]=V t 

, „ , , . , 315VTTv 1155VT7v 3 9009VTTV 5 , — 7 r — 

Define the interval (note that for Legendre polynomials : v x — — 71 7 + 7 25740V 19 v' + 48620 V 19 v 

l;v 2 =l) 65 


64 


vl=- 1 / 2 ;v 2 = 1 / 2 ; 
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Check Inner Product: 


MatrixForm[Table[ip[i,j], {i,0,terms-l}, b Aterms-1}]] 

5 


( 1 

0 

0 

0 

0 

0 

0 

0 

0 

0> 




0 

1 

0 

0 

0 

0 

0 

0 

0 

0 




0 

0 

1 

0 

0 

0 

0 

0 

0 

0 




0 

0 

0 

1 

0 

0 

0 

0 

0 

0 


10 


0 

0 

0 

0 

1 

0 

0 

0 

0 

0 



0 

0 

0 

0 

0 

1 

0 

0 

0 

0 




0 

0 

0 

0 

0 

0 

1 

0 

0 

0 




0 

0 

0 

0 

0 

0 

0 

1 

0 

0 




0 

0 

0 

0 

0 

0 

0 

0 

1 

0 
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k 0 

0 

0 

0 

0 

0 

0 

0 

0 

1 , 




We claim: 

1. A method of sampling a band-limited function, the 20 
method comprising: 

generating transformed data using a discrete Fourier trans- 
form which samples a band-limited function f(x) with- 
out interpolating or modifying received data associated 
with the function f(x), wherein an interval between 25 
repeated copies in a periodic extension of the function 
f(x) obtained from the discrete Fourier transform is asso- 
ciated with a sampling ratio Q, defined as a ratio of a 
sampling frequency to a band-limited frequency, and 
wherein Q is assigned a value between 1 and 2 such that 30 
substantially no aliasing occurs in the transformed data; 
and 

retrieving a phase in the received data based on the trans- 
formed data. 

2. The method of claim 1 wherein the band-limited func- 35 
tion is a sinc(x) function. 

3. The method of claim 1, wherein Q=(lambda*f/#)/dx and 
lambda is a monochromatic wavelength of a scalar electro- 
magnetic field, f/# is an f-number of an optical system and dx 

is an image plane sampling interval. 40 

4. The method of claim 1, wherein the received data is from 
a confined circular aperture of diameter D. 

5. The method of claim 1, further comprising adjusting the 
optical system based on the feedback. 

6 . The method of claim 1, wherein the phase is used as 45 
feedback to an optical system. 

7. A system for sampling a band-limited function, the sys- 
tem comprising: 

a processor; 

a first processor controlling the processor to generate trans- 50 
formed data using a discrete Fourier transform which 
samples a band-limited function f(x) without interpolat- 
ing or modifying received data associated with the func- 
tion f(x), wherein an interval between repeated copies in 
a periodic extension of the function f(x) obtained from 
the discrete Fourier transform is associated with a sam- 


18 

pling ratio Q, defined as a ratio of a sampling frequency 
to a band-limited frequency, and wherein Q is assigned a 
value between 1 and 2 such that substantially no aliasing 
occurs in the transformed data; and 

a second module controlling the processor to retrieve a 
phase in the received data based on the transformed data. 

8. The system of claim 7, wherein the band-limited func- 
tion is a sinc(x) function. 

9. The system of claim 7, wherein Q=(l ambda * f/ # )/ dx and 
lambda is a monochromatic wavelength of a scalar electro- 
magnetic field, f/# is an f-number of an optical system and dx 
is an image plane sampling interval. 

10. The system of claim 7, wherein the received data is 
from a confined circular aperture of diameter D. 

11. The system of claim 7, further comprising a third mod- 
ule controlling the processor to adjust the optical system 
based on the feedback. 

12. The system of claim 7, wherein the phase is used as 
feedback to an optical system. 

13. A non-transitory computer-readable storage medium 
storing instructions which, when executed by a computing 
device, cause the computing device to sample a band-limited 
function, the instructions comprising: 

generating transformed data using a discrete Fourier trans- 
form which samples a band- limited function f(x) with- 
out interpolating or modifying received data associated 
with the function f(x), wherein an interval between 
repeated copies in a periodic extension of the function 
f(x) obtained from the discrete Fourier transform is asso- 
ciated with a sampling ratio Q, defined as a ratio of a 
sampling frequency to a band-limited frequency, and 
wherein Q is assigned a value between 1 and 2 such that 
substantially no aliasing occurs in the transformed data; 
and 

retrieving a phase in the received data based on the trans- 
formed data. 

14. The non-transitory computer-readable storage medium 
of claim 13, wherein the band-limited function is a sinc(x) 
function. 

15. The non-transitory computer-readable storage medium 
of claim 13, wherein Q=(lambda * f/ #)/ dx and lambda is a 
monochromatic wavelength of a scalar electromagnetic field, 
f/# is an f-number of an optical system and dx is an image 
plane sampling interval. 

16. The non-transitory computer-readable storage medium 
of claim 13, wherein the received data is from a confined 
circular aperture of diameter D. 

17. The non-transitory computer-readable storage medium 
of claim 13, the instructions further comprising adjusting the 
optical system based on the feedback. 

18. The non-transitory computer-readable storage medium 
of claim 13, wherein the phase is used as feedback to an 
optical system. 





